#############################################################################
##
#W  refine.gi                DifSets Package                     Dylan Peifer
##
##  Functions find all possible difference set/sum preimages of difference
##  sums.
##

#############################################################################
##
#F  IsCoPrime(<a>, <b>)
##
#  判断给定的两个元素a和b是否互素
InstallGlobalFunction( IsCoPrime, function (a, b)
    if Gcd(a, b)=1 then
        return true;
    else
        return false;
    fi;
end );

#############################################################################
##
#F  AllCoPrime(<v>)
##
#  给出小于v的数中与v互素的所有数
InstallGlobalFunction( AllCoPrime, function (v)
    local coprimes, i;
    coprimes := [];
    for i in [2..v-1] do
        if IsCoPrime(v, i) then
            Add(coprimes, i);
        fi;
    od;
    return coprimes;
end );


#############################################################################
##
#F  Multipliers
##
InstallGlobalFunction( Multipliers, function (G, k, lambda)
    local allcoprimes, n, p, possibleP, flag;
    allcoprimes := AllCoPrime(Size(G));  # 2..v-1中所有与v互素的数

    # 判断是否为n=k-lambda的素因子
    possibleP := [];
    n := k-lambda;
    flag := false;
    for p in allcoprimes do
        if IsPrime(p) and n mod p =0 then
            Add(possibleP, p);
            if p > lambda then
                flag := true;
            fi;
        fi;
    od;

    # p > lambda
    if flag then
        for p in possibleP do
            if not p > lambda then
                Filtered(possibleP, x->(x<>p));
            fi;
        od;
    fi;

    # 即不满足乘子定理 也不满足乘子猜想 直接返回所有与v互素的数
    if IsEmpty(possibleP) then
        # Print("使用所有与v互素的数作为乘子", "\n");
        return allcoprimes;
    fi;
    return possibleP;
end );

#############################################################################
##
#F  MultiplierTheorem
##
InstallGlobalFunction( MultiplierTheorem, function (G, k, lambda)
    local allcoprimes, n, p, possibleP, flag;
    allcoprimes := AllCoPrime(Size(G));  # 2..v-1中所有与v互素的数

    # 判断是否为n=k-lambda的素因子
    possibleP := [];
    n := k-lambda;
    flag := false;
    for p in allcoprimes do
        if IsPrime(p) and n mod p =0 then
            Add(possibleP, p);
            if p > lambda then
                flag := true;
            fi;
        fi;
    od;

    # p > lambda
    if flag then
        # 满足乘子定理
        for p in possibleP do
            if not p > lambda then
                Filtered(possibleP, x->(x<>p));
            fi;
        od;
        return possibleP;
    else
        # 仅满足乘子猜想，返回空列表
        return [];
    fi;
end );

#############################################################################
##
#F  Export
##
InstallGlobalFunction( Export, function ()
    local order, orders;
    orders:=[];
    for order in [2..100] do
        if(Size(FindAllMultipliersOfDifferenceSet(CyclicGroup(order)))>0) then
            Add(orders, order);
        fi;
    od;
    return orders;
end );

#############################################################################
##
#F  FindAllMultipliersOfDifferenceSet
##
InstallGlobalFunction( FindAllMultipliersOfDifferenceSet, function (G)
    local v, k, lambda, multipliers, part;
    v := Size(G);
    multipliers := [];
    for k in PossibleDifferenceSetSizes(G) do
        lambda := (k*k - k) / (v-1);
        part := MultiplierTheorem(G, k, lambda);
        Append(multipliers, part);
    od;

    return multipliers;
end );


#############################################################################
##
#F  TracksOfGroup( <G>, <multiplier> )
##
#  将群划分为多个轨道
InstallGlobalFunction( TracksOfGroup, function (G, multiplier)
    local tracksOfGroup, v, elements, usedmap, next, newtrack, ele, pos;
    tracksOfGroup := [];  # will store final results.
    v:=Size(G);
    elements := Elements(G);
    usedmap := List([1..v], x->1);  # 1 表示暂未划分到轨道中
    next := 1;  # next 表示下一个待划分的元素的索引
    while next <= v do
        # 开始一个新的轨道
        newtrack := [];
        ele := elements[next];
        usedmap[next] := 0;
        pos := next;
        while not pos in newtrack do
            Add(newtrack, pos);
            usedmap[pos] := 0;
            ele := ele^multiplier;
            pos := Position(elements, ele);
        od;

        # 当前轨道结束
        Add(tracksOfGroup, newtrack);

        while next <= v and usedmap[next] = 0 do
            next := next + 1;
        od;

    od;
    return tracksOfGroup;
end );

#############################################################################
###
##F  TracksOfCoset( <cosets>, <tracksOfGroup> )
###
InstallGlobalFunction( TracksOfCoset, function (cosets, tracksOfGroup)
    local i, coset, track, tracksOfCoset, newtrack;
    tracksOfCoset := List([1..Size(cosets)], x->[]);  # will store the final result
    for i in [1..Size(cosets)] do
        coset := cosets[i];
        for track in tracksOfGroup do
            # 计算交集
            newtrack := Intersection(coset, track);
            if not IsEmpty(newtrack) then
                Add(tracksOfCoset[i], newtrack);
            fi;
        od;
    od;

    return tracksOfCoset;
end );


#############################################################################
##
#F  AllRefineDifferenceSets( <G>, <N>, <difsums> )
##
InstallGlobalFunction( AllRefinedDifferenceSets, function (G, N, difsums)
    local v, k, lambda, dif, cosets, difsets, S, opts, opt, D;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);

    difsets := []; # will store found difsets

    for S in difsums do

        # iterate through preimages D of S
        opts := List([1..Size(cosets)], x->Combinations(cosets[x], S[x]));
        for opt in IteratorOfCartesianProduct(opts) do
            D := Flat(opt);

            # if D is a difset then add it to list
            if IsDifferenceSetByTable(dif, D, v, lambda) then
                Add(difsets, D);
            fi;

        od;

    od;

    return difsets;
end );

#############################################################################
##
#F  NrAllRefinedSets( <G>, <N>, <difsums> )
##
InstallGlobalFunction( NrAllRefinedSets, function (G, N, difsums)
    return Sum(difsums, S->Product(S, x->Binomial(Size(N), x)));
end );

#############################################################################
##
#F  SomeRefinedDifferenceSets( <G>, <N>, <difsums> )
##
InstallGlobalFunction( SomeRefinedDifferenceSets, function (G, N, difsums)
    local v, k, lambda, dif, cosets, difsets, S, opts, opt, D;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);

    difsets := []; # will store found difsets

    for S in difsums do

        # iterate through preimages D of S, forcing choice of identity
        opts := List([1..Size(cosets)], x->Combinations(cosets[x], S[x]));
        if S[1] > 0 then opts[1] := Filtered(opts[1], x->(1 in x)); fi;

        for opt in IteratorOfCartesianProduct(opts) do
            D := Flat(opt);

            # if D is a difset then add it to list
            if IsDifferenceSetByTable(dif, D, v, lambda) then
                Add(difsets, D);
            fi;

        od;

    od;

    return difsets;
end );

#############################################################################
##
#F  SubsetSum( <arr>, <sum> )
##
InstallGlobalFunction( SubsetSum, function (nums, k)
    local res, subset, findSubsets, i;
    res := [];
    subset := [];
    findSubsets := function(nums, k, res, subset, index)
        if k = 0 then
            Add(res, ShallowCopy(subset));
            return;
        fi;
        for i in [index..Size(nums)] do
            Add(subset, i);
            findSubsets(nums, k - nums[i], res, subset, i + 1);
            Remove(subset);
        od;
    end;
    findSubsets(nums, k, res, subset, 1);
    return res;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSetsWithMultiplier( <G>, <N>, <difsums> )
##
InstallGlobalFunction( SomeRefinedDifferenceSetsWithMultiplier, function (G, N, difsums)
    local v, k, lambda, dif, cosets, difsets, S, opts, opt, D, multipliers, multiplier, tracksOfGroup, tracksOfCoset;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);
    # Print("SomeRefinedDifferenceSetsWithMultiplier>>>cosets=",cosets, "\n");

    difsets := []; # will store found difsets

    # if IsCyclic(G) and Size(multipliers) > 0 then
    if IsAbelian(G) then
        multipliers := Multipliers(G, k, lambda);
        # 使用轨道
        # Print("SomeRefinedDifferenceSetsWithMultiplier>>>使用轨道","\n");
        for multiplier in multipliers do
            tracksOfGroup := TracksOfGroup(G, multiplier);
            # Print("SomeRefinedDifferenceSetsWithMultiplier>>>tracksOfGroup=",tracksOfGroup, "\n");
            tracksOfCoset := TracksOfCoset(cosets, tracksOfGroup);
            # Print("SomeRefinedDifferenceSetsWithMultiplier>>>tracksOfCoset=",tracksOfCoset, "\n");
            for S in difsums do
                # iterate through preimages D of S
                # Print("SomeRefinedDifferenceSetsWithMultiplier>>>S=",S, "\n");
                opts := List([1..Size(cosets)], x->SubsetSum(List([1..Size(tracksOfCoset[x])], y-> Size(tracksOfCoset[x][y])), S[x]));
                # Print("SomeRefinedDifferenceSetsWithMultiplier>>>索引opts=",opts, "\n");
                opts := List([1..Size(cosets)], x->List(opts[x], y->List(y, z-> tracksOfCoset[x][z])));
                # Print("SomeRefinedDifferenceSetsWithMultiplier>>>元素opts=",opts, "\n");
                opts:= List(opts, x->List(x, y->Flat(y)));
                # Print("SomeRefinedDifferenceSetsWithMultiplier>>>展平opts=",opts, "\n");
                # FIXME 需过滤！！！
                if S[1] > 0 then opts[1] := Filtered(opts[1], x->(1 in x)); fi;
                # Print("SomeRefinedDifferenceSetsWithMultiplier>>>过滤opts=",opts, "\n");
                for opt in IteratorOfCartesianProduct(opts) do
                    D := Flat(opt);
                    # if D is a difset then add it to list
                    if IsDifferenceSetByTable(dif, D, v, lambda) then
                        Add(difsets, D);
                    fi;
                od;
            od;
        od;
    else
        # as usual
        for S in difsums do
            # iterate through preimages D of S, forcing choice of identity
            opts := List([1..Size(cosets)], x->Combinations(cosets[x], S[x]));
            if S[1] > 0 then opts[1] := Filtered(opts[1], x->(1 in x)); fi;
            for opt in IteratorOfCartesianProduct(opts) do
                D := Flat(opt);
                # if D is a difset then add it to list
                if IsDifferenceSetByTable(dif, D, v, lambda) then
                    Add(difsets, D);
                fi;
            od;
        od;
    fi;
    return difsets;
end );


#############################################################################
##
#F  SomeRefinedDifferenceSetsUseMultiplier( <G>, <N>, <difsums> )
##
InstallGlobalFunction( SomeRefinedDifferenceSetsUseMultiplier, function (G, N, difsums, k)
    local v, lambda, Q, dif, cosets, cosetsInTracks, T, difsum, opts,opt, difset, difsets, multiplier;
    v:= Size(G);
    lambda := k*(k-1)/(v-1);
    Q:=G/N;
    dif := DifferenceTable(G);
    cosets:=CosetsTable(G,N);
    multiplier := Multipliers(G,k,lambda)[1];
    T:=TracksOfGroup(Q,multiplier);
    cosetsInTracks := List(T, x->List(x,y->cosets[y]));
    cosetsInTracks := List(cosetsInTracks, x->Flat(x));
    difsets:=[];
    for difsum in difsums do
        opts := List([1..Size(difsum)], x->Combinations(cosetsInTracks[x], difsum[x]));
        for opt in IteratorOfCartesianProduct(opts) do
            difset := Flat(opt);
            if IsDifferenceSetByTable(dif, difset, v, lambda) then
                Add(difsets, difset);
            fi; 
        od;
    od;
    return difsets;
end );

#############################################################################
##
#F  NrSomeRefinedSets( <G>, <N>, <difsums> )
##
InstallGlobalFunction( NrSomeRefinedSets, function (G, N, difsums)
    local f;

    f := function(S)
        if S[1] > 0 then
            return Binomial(Size(N)-1, S[1]-1)
                   * Product(S{[2..Length(S)]}, x->Binomial(Size(N), x));
        else
            return Product(S, x->Binomial(Size(N), x));
        fi;
    end;

    return Sum(difsums, S->f(S));
end );

#############################################################################
##
#F  SomeRefinedDifferenceSetPairs( <G>, <N>, <difsums>, <e>, <lambda>)
##
# FEATURE SomeRefinedDifferenceSetPairs(G, N, difsums, e, lambda)
InstallGlobalFunction( SomeRefinedDifferenceSetPairs, function (G, N, difsums, e, lambda)
    local v, dif, cosets, difsetpairs, S, opts1, opts2, opt1, opt2, D1, D2;
    if Length(difsums) = 0 then return []; fi;
    v := Size(G);
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);
    difsetpairs := []; # will store found difsetpairs

    # Print("SomeRefinedDifferenceSetPairs>>>|G|=", v, ", |N|=", Size(N), "difsums=", difsums, "\n");
    # Print("SomeRefinedDifferenceSetPairs>>>cosets=", cosets, "\n");
    for S in difsums do
        # iterate through preimages D of S, forcing choice of identity
        opts1 := List([1..Size(cosets)], x->Combinations(cosets[x], S[1][x]));
        opts2 := List([1..Size(cosets)], x->Combinations(cosets[x], S[2][x]));
        if e > 0 then  # forcing choice of identity when e > 0
            opts1[1] := Filtered(opts1[1], x->(1 in x));
            opts2[1] := Filtered(opts2[1], x->(1 in x));
        fi;

        # if e > 1 then  # forcing choice of two exact element when e > 1
        #     # giving another element to filter out extra options, at this point, the identity element must exist.
        #     opts1[1] := Filtered(opts1[1], x->(2 in x));
        #     opts2[1] := Filtered(opts2[1], x->(2 in x));
        # fi;

        for opt1 in IteratorOfCartesianProduct(opts1) do
            D1 := Flat(opt1);
            for opt2 in IteratorOfCartesianProduct(opts2) do
                D2 := Flat(opt2);
                # check if the length of intersection equals to e
                if Length(Intersection(D1,D2)) <> e then
                    continue;
                fi;
                # if D1 D2 is a difsetpair then add it to list
                if IsDifferenceSetPairByTable(dif, D1, D2, v, e, lambda) then
                    Add(difsetpairs, [D1,D2]);
                    break; # break when we found a dsp
                fi;
            od;
        od;

    od;
    # Print("SomeRefinedDifferenceSetPairs>>>difsetpairs=", difsetpairs, "\n");
    return difsetpairs;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSetPairsWithMultiplier( <G>, <N>, <difsums>, <e>, <lambda>)
##
# FEATURE SomeRefinedDifferenceSetPairsWithMultiplier(G, N, difsums, e, lambda)
InstallGlobalFunction( SomeRefinedDifferenceSetPairsWithMultiplier, function (G, N, difsums, e, lambda)
    local v, dif, cosets, difsetpairs, S, opts1, opts2, opt1, opt2, D1, D2, multiplier, multipliers, tracksOfGroup, tracksOfCoset;
    if Length(difsums) = 0 then return []; fi;
    v := Size(G);
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);
    difsetpairs := []; # will store found difsetpairs

    # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>|G|=", v, ", |N|=", Size(N), "difsums=", difsums, "\n");
    # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>cosets=", cosets, "\n");
    if IsAbelian(G) then
        multipliers := Multipliers(G, e, lambda);
        # 使用轨道
        # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>使用轨道","\n");
        for multiplier in multipliers do
            tracksOfGroup := TracksOfGroup(G, multiplier);
            # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>tracksOfGroup=",tracksOfGroup, "\n");
            tracksOfCoset := TracksOfCoset(cosets, tracksOfGroup);
            # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>tracksOfCoset=",tracksOfCoset, "\n");

            
            for S in difsums do
                # iterate through preimages D of S
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>S=",S, "\n");
                opts1 := List([1..Size(cosets)], x->SubsetSum(List([1..Size(tracksOfCoset[x])], y-> Size(tracksOfCoset[x][y])), S[1][x]));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>索引opts1=",opts1, "\n");
                opts1 := List([1..Size(cosets)], x->List(opts1[x], y->List(y, z-> tracksOfCoset[x][z])));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>元素opts1=",opts1, "\n");
                opts1:= List(opts1, x->List(x, y->Flat(y)));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>展平opts1=",opts1, "\n");
    
                opts2 := List([1..Size(cosets)], x->SubsetSum(List([1..Size(tracksOfCoset[x])], y-> Size(tracksOfCoset[x][y])), S[2][x]));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>索引opts2=",opts2, "\n");
                opts2 := List([1..Size(cosets)], x->List(opts2[x], y->List(y, z-> tracksOfCoset[x][z])));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>元素opts2=",opts2, "\n");
                opts2:= List(opts2, x->List(x, y->Flat(y)));
                # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>展平opts2=",opts2, "\n");

                if e > 0 then  # forcing choice of identity when e > 0
                    opts1[1] := Filtered(opts1[1], x->(1 in x));
                    opts2[1] := Filtered(opts2[1], x->(1 in x));
                fi;

    
                for opt1 in IteratorOfCartesianProduct(opts1) do
                    D1 := Flat(opt1);
                    for opt2 in IteratorOfCartesianProduct(opts2) do
                        D2 := Flat(opt2);
                        # check if the length of intersection equals to e
                        if Length(Intersection(D1,D2)) <> e then
                            continue;
                        fi;
                        # if D1 D2 is a difsetpair then add it to list
                        if IsDifferenceSetPairByTable(dif, D1, D2, v, e, lambda) then
                            Add(difsetpairs, [D1,D2]);
                            break; # break when we found a dsp
                        fi;
                    od;
                od;
        
            od;
            
        od;

        
    else
        for S in difsums do
            # iterate through preimages D of S, forcing choice of identity
            opts1 := List([1..Size(cosets)], x->Combinations(cosets[x], S[1][x]));
            opts2 := List([1..Size(cosets)], x->Combinations(cosets[x], S[2][x]));
            if e > 0 then  # forcing choice of identity when e > 0
                opts1[1] := Filtered(opts1[1], x->(1 in x));
                opts2[1] := Filtered(opts2[1], x->(1 in x));
            fi;
    
            # if e > 1 then  # forcing choice of two exact element when e > 1
            #     # giving another element to filter out extra options, at this point, the identity element must exist.
            #     opts1[1] := Filtered(opts1[1], x->(2 in x));
            #     opts2[1] := Filtered(opts2[1], x->(2 in x));
            # fi;
    
            for opt1 in IteratorOfCartesianProduct(opts1) do
                D1 := Flat(opt1);
                for opt2 in IteratorOfCartesianProduct(opts2) do
                    D2 := Flat(opt2);
                    # check if the length of intersection equals to e
                    if Length(Intersection(D1,D2)) <> e then
                        continue;
                    fi;
                    # if D1 D2 is a difsetpair then add it to list
                    if IsDifferenceSetPairByTable(dif, D1, D2, v, e, lambda) then
                        Add(difsetpairs, [D1,D2]);
                        break; # break when we found a dsp
                    fi;
                od;
            od;
    
        od;
    fi;
    # Print("SomeRefinedDifferenceSetPairsWithMultiplier>>>difsetpairs=", difsetpairs, "\n");
    return difsetpairs;
end );


#############################################################################
##
#F  SomeRefinedDivisibleDifferenceSets( <G>, <N>, <difsums>, <w>, <v>, <l>, <lambda1>, <lambda2> )
##
# FEATURE SomeRefinedDivisibleDifferenceSets(G, N, difsums, w, v, l, lambda1, lambda2)
InstallGlobalFunction( SomeRefinedDivisibleDifferenceSets, function (G, N, difsums, w, v, l, lambda1, lambda2)
    local dif, cosets, dividifsets, S, opts, opt, D;
    if Length(difsums) = 0 then return []; fi;
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);

    dividifsets := []; # will store found dividifsets

    for S in difsums do
        # iterate through preimages D of S, forcing choice of identity
        opts := List([1..Size(cosets)], x->Combinations(cosets[x], S[x]));
        if S[1] > 0 then opts[1] := Filtered(opts[1], x->(1 in x)); fi;
        for opt in IteratorOfCartesianProduct(opts) do
            D := Flat(opt);
            if IsDivisibleDifferenceSetByTable(G, dif, D, w, v, l, lambda1, lambda2) then
                Add(dividifsets, D);
            fi;

        od;

    od;
    return dividifsets;
end );


#############################################################################
##
#F  SomeRefinedDivisibleDifferenceSetPairs( <G>, <N>, <difsums>, <w>, <v>, <l1>, <l2>, <le>, <lambda1>, <lambda2> )
##
# FEATURE SomeRefinedDivisibleDifferenceSetPairs(G, N, difsums, w, v, l1, l2, le, lambda1, lambda2)
InstallGlobalFunction( SomeRefinedDivisibleDifferenceSetPairs, function (G, N, difsums, w, v, l1, l2, le, lambda1, lambda2)
    local dif, cosets, dividifsetpairs, S, opts1, opts2, opt1, opt2, D1, D2;
    if Length(difsums) = 0 then return []; fi;
    dif := DifferenceTable(G);
    cosets := CosetsTable(G, N);
    dividifsetpairs := []; # will store found dividifsetpairs
    for S in difsums do
        # iterate through preimages D of S, forcing choice of identity
        opts1 := List([1..Size(cosets)], x->Combinations(cosets[x], S[1][x]));
        opts2 := List([1..Size(cosets)], x->Combinations(cosets[x], S[2][x]));
        for opt1 in IteratorOfCartesianProduct(opts1) do
            D1 := Flat(opt1); 
            for opt2 in IteratorOfCartesianProduct(opts2) do
                D2 := Flat(opt2);
                if IsDivisibleDifferenceSetPairByTable(G, dif, D1, D2, w, v, l1, l2, le, lambda1, lambda2) then
                    Add(dividifsetpairs, [D1,D2]);
                fi;
            od;
        od;

    od;

    return dividifsetpairs;
end );

#############################################################################
##
#F  IsDifferenceSetByTable( <table>, <D>, <v>, <lambda> )
##
InstallGlobalFunction( IsDifferenceSetByTable, function (table, D, v, lambda)
    # method is identical to IsDifferenceSet but uses precomputed values
    local count, i, j, index;

    count := List([1..v], x->0);

    for i in D do
    for j in D do

        index := table[i][j];
        count[index] := count[index] + 1;

        if not index = 1 and count[index] > lambda then
            return false;
        fi;

    od;
    od;

    return true;
end );

#############################################################################
##
#F  IsDifferenceSetPairByTable( <table>, <D1>, <D2>, <v>, <e>, <lambda> )
##
# FEATURE IsDifferenceSetPairByTable(table, D1, D2, v, e, lambda)
InstallGlobalFunction( IsDifferenceSetPairByTable, function (table, D1, D2, v, e, lambda)
    local diff, i, j, index;
    diff := List([1..v], x->0);

    for i in D1 do
        for j in D2 do
            index := table[i][j];
            diff[index] := diff[index] + 1;

            if not index = 1 and diff[index] > lambda then
                return false;
            fi;
        od;
    od;
    if diff[1] <> e then
        return false;
    fi;
    return true;
end );

#############################################################################
##
#F  IsDivisibleDifferenceSetByTable(<G>, <table>, <D, <w>, <v>, <l>, <lambda1>, <lambda2>)
##
# FEATURE IsDivisibleDifferenceSetByTable(G, table, D, w, v, l, lambda1, lambda2)
InstallGlobalFunction( IsDivisibleDifferenceSetByTable, function (G, table, D, w, v, l, lambda1, lambda2)
    local n, diff, i, j, index, Hs, H, eles;
    Hs := NormalSubgroups(G);
    for i in Hs do
        if Size(i)=v then
            H := i;
            break;
        fi;
    od;

    n := w*v;
    diff := List([1..n], x->0);
    for i in D do
        for j in D do
            index := table[i][j];
            diff[index] := diff[index] +1 ;
        od;
    od;

    eles := Elements(G);
    for i in [2..Length(eles)] do # ignore identity element of G
        i := eles[i];
        if i in H then
            if  diff[Position(eles, i)] <> lambda1 then
                return false;
            fi;
        else
            if  diff[Position(eles, i)] <> lambda2 then
                return false;
            fi;
        fi;
    od;

    return true;
end );

#############################################################################
##
#F  IsDivisibleDifferenceSetPairByTable(<G>, <table>, <D1>, <D2>, <w>, <v>, <l1>, <l2>, <le>, <lambda1>, <lambda2>)
##
# FEATURE IsDivisibleDifferenceSetPairByTable(G, table, D1, D2, w, v, l1, l2, le, lambda1, lambda2)
InstallGlobalFunction( IsDivisibleDifferenceSetPairByTable, function (G, table, D1, D2, w, v, l1, l2, le, lambda1, lambda2)
    local n, diff, i, j, index, Hs, H, eles;
    Hs := NormalSubgroups(G);
    for i in Hs do
        if Size(i)=v then
            H := i;
            break;
        fi;
    od;

    n := w*v;
    diff := List([1..n], x->0);
    for i in D1 do
        for j in D2 do
            index := table[i][j];
            diff[index] := diff[index] +1 ;
        od;
    od;

    eles := Elements(G);
    for i in [2..Length(eles)] do # ignore identity element of G
        i := eles[i];
        if i in H then
            if  diff[Position(eles, i)] <> lambda1 then
                return false;
            fi;
        else
            if  diff[Position(eles, i)] <> lambda2 then
                return false;
            fi;
        fi;
    od; 
    return true;
end );

#############################################################################
##
#F  AllRefinedDifferenceSums( <G>, <N1>, <N2>, <difsums> )
##
InstallGlobalFunction( AllRefinedDifferenceSums, function (G, N1, N2, difsums)
    local v, k, lambda, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts, opt, perms, perm, s;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);

    # compute some numbers for function calls below
    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    id := lambda*Size(N2) + k - lambda;
    nonid := lambda*Size(N2);

    newdifsums := []; # will store found difsums

    for S in difsums do

        # iterate through all partitions of S into new cosets
        opts := DifferenceSumPreImagesOptions(u3, w2,  S);
        for opt in IteratorOfCartesianProduct(opts) do

            # testing identity coeff can be done before permutations
            if not Sum(Flat(opt), x->x^2) = id then
                continue;
            fi;

            # iterate through all permutations of this partition
            perms := DifferenceSumPreImagesPermutations(opt);
            for perm in IteratorOfCartesianProduct(perms) do

                # if difference sum is found then add to list
                s := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm);
                if IsDifferenceSumByTable(dif, s, u2, nonid) then
                    Add(newdifsums, s);
                fi;

            od;

        od;

    od;

    return newdifsums;
end );

#############################################################################
##
#F  NrAllRefinedSums( <G>, <N1>, <N2>, <difsums> )
##
InstallGlobalFunction( NrAllRefinedSums, function (G, N1, N2, difsums)
    local v, k, lambda, total, S, opts, opt, perms;

    if Length(difsums) = 0 then return 0; fi;

    v := Size(G);
    k := Sum(difsums[1]);
    lambda := k*(k-1)/(v-1);

    total := 0;
    for S in difsums do
        opts := DifferenceSumPreImagesOptions(Size(N1/N2), Size(N2), S);

        for opt in IteratorOfCartesianProduct(opts) do
            if not Sum(Flat(opt), x->x^2) = lambda*Size(N2)+k-lambda then
                continue;
            fi;

            perms := DifferenceSumPreImagesPermutations(opt);
            total := total + Product(perms, x->Length(x));

        od;

    od;

    return total;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSums( <G>, <N1>, <N2>, <difsums> )
##
InstallGlobalFunction( SomeRefinedDifferenceSums, function (G, N1, N2, difsums)
    local v, k, lambda, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts, opt, perms, perm, s;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);

    # compute some numbers for function calls below
    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    id := lambda*Size(N2) + k - lambda;
    nonid := lambda*Size(N2);

    newdifsums := []; # will store found difsums

    for S in difsums do

        # iterate through all partitions of S into new cosets
        opts := DifferenceSumPreImagesOptions(u3, w2,  S);
        for opt in IteratorOfCartesianProduct(opts) do

            # testing identity coeff can be done before permutations
            if not Sum(Flat(opt), x->x^2) = id then
                # Print("过滤掉组合（未全排列）【",opt,"】\n");
                continue;
            fi;

            # iterate through all permutations of this partition
            perms := DifferenceSumPreImagesPermutationsForced(opt);
            for perm in IteratorOfCartesianProduct(perms) do

                # if difference sum is found then add to list
                s := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm);
                if IsDifferenceSumByTable(dif, s, u2, nonid) then
                    Add(newdifsums, s);
                else
                    # Print("过滤掉差集和【",s,"】\n");
                fi;

            od;

        od;

    od;

    return newdifsums;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSumsOptimized( <G>, <N1>, <N2>, <difsums> )
##
InstallGlobalFunction( SomeRefinedDifferenceSumsOptimized, function (G, N1, N2, difsums)
    local v, k, lambda, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts, opt, perms, perm, s, multiplier, tracksOfGroup, track, flag, parts;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    k := Sum(difsums[1]); # assume all difsums have same k
    lambda := k*(k-1)/(v-1);
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);

    # compute some numbers for function calls below
    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    id := lambda*Size(N2) + k - lambda;
    nonid := lambda*Size(N2);

    newdifsums := []; # will store found difsums

    # optimization：在同一轨道的陪集，其差集和元素相同
    multiplier := Multipliers(G,k,lambda)[1];
    tracksOfGroup := TracksOfGroup(G/N2, multiplier);

    for S in difsums do

        # iterate through all partitions of S into new cosets
        opts := DifferenceSumPreImagesOptions(u3, w2,  S);
        for opt in IteratorOfCartesianProduct(opts) do

            # testing identity coeff can be done before permutations
            if not Sum(Flat(opt), x->x^2) = id then
                # Print("过滤掉组合（未全排列）【",opt,"】\n");
                continue;
            fi;

            # iterate through all permutations of this partition
            perms := DifferenceSumPreImagesPermutationsForced(opt);
            for perm in IteratorOfCartesianProduct(perms) do

                # if difference sum is found then add to list
                s := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm);
                flag := true;  # 该差集和满足在同一轨道下，元素相同
                for track in tracksOfGroup do
                    if Length(track)=1 then
                        continue;
                    else
                        parts := s{track};  # 该轨道下对应的差集和元素，理论上应该相等
                        if(Length(Set(parts))>1) then
                            flag := false;
                            break;
                        fi;
                    fi;
                od;
                if flag = true and IsDifferenceSumByTable(dif, s, u2, nonid) then
                    Add(newdifsums, s);
                fi;

            od;

        od;

    od;
    return newdifsums;
end );

#############################################################################
##
#F  TracksToTracksTable(<G>, <N1>, <N2>)
##
InstallGlobalFunction( TracksToTracksTable, function (G, N1, N2, multiplier)
    local Q1, Q2, Tracks1, Tracks2, Cosets1, Cosets2, results, idxi, idxj, compo, curTrack;
    Q1 := G/N1;  # 商群
    Q2 := G/N2;  # 商群

    Tracks1 := TracksOfGroup(Q1, multiplier);
    Tracks2 := TracksOfGroup(Q2, multiplier);

    Cosets1 := CosetsTable(G,N1);
    Cosets2 := CosetsTable(G,N2);

    Tracks1:=List(Tracks1, x->List(x, y->Cosets1[y]));
    Tracks2:=List(Tracks2, x->List(x, y->Cosets2[y]));

    Tracks1:=List(Tracks1, x->Flat(x));
    Tracks2:=List(Tracks2, x->Flat(x));

    results:=[];
    for idxi in [1..Size(Tracks1)] do
        curTrack := Tracks1[idxi];
        compo := [];
        for idxj in [1..Size(Tracks2)] do
            if(Tracks2[idxj][1] in curTrack) then
                Add(compo, idxj);
            fi;
        od;
        Add(results, compo);
    od;

    return results;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSumsUseMultiplier( G, k, lambda, N1, N2, difsums )
##
InstallGlobalFunction( SomeRefinedDifferenceSumsUseMultiplier, function (G, k, lambda, N1, N2, difsums)
    local v, tracksMap, newdifsums, multiplier, S, opts, opt, perms, perm, idxi, idxj, s;

    # handle special case of empty input
    if Length(difsums) = 0 then return []; fi;

    v := Size(G);
    # Print("v=",v,"\n");
    multiplier := Multipliers(G,k,lambda)[1];
    # Print("multiplier=",multiplier,"\n");
    tracksMap := TracksToTracksTable(G, N1, N2,multiplier);
    # Print("tracksMap=",tracksMap,"\n");

    newdifsums := []; # will store found difsums

    for S in difsums do
        opts :=List([1..Size(S)], x->RestrictedPartitions(S[x]+Size(tracksMap[x]),[1..S[x]+Size(tracksMap[x])], Size(tracksMap[x])));
        # Print("opts=",opts,"\n");
        opts:=List(opts, x->List(x,y->List(y, z->z-1)));
        # Print("opts=",opts,"\n");
        
        for opt in IteratorOfCartesianProduct(opts) do
            # Print("opt=",opt,"\n");
            perms:= List(opt,x->PermutationsList(x));
            # Print("perms=",perms,"\n");
            for perm in IteratorOfCartesianProduct(perms) do
                # translate to corresponding position
                # Print("perm=",perm,"\n");
                s:=[];
                for idxi in [1..Size(tracksMap)] do
                    for idxj in [1..Size(tracksMap[idxi])] do
                        s[tracksMap[idxi][idxj]]:=perm[idxi][idxj];
                    od;
                od;

                Add(newdifsums, s);
            od;
        od;
    od;

    return newdifsums;
end );

#############################################################################
##
#F  NrSomeRefinedSums( <G>, <N1>, <N2>, <difsums> )
##
InstallGlobalFunction( NrSomeRefinedSums, function (G, N1, N2, difsums)
    local v, k, lambda, total, S, opts, opt, perms;

    if Length(difsums) = 0 then return 0; fi;

    v := Size(G);
    k := Sum(difsums[1]);
    lambda := k*(k-1)/(v-1);

    total := 0;
    for S in difsums do
        opts := DifferenceSumPreImagesOptions(Size(N1/N2), Size(N2), S);

        for opt in IteratorOfCartesianProduct(opts) do
            if not Sum(Flat(opt), x->x^2) = lambda*Size(N2)+k-lambda then
                continue;
            fi;

            perms := DifferenceSumPreImagesPermutationsForced(opt);
            total := total + Product(perms, x->Length(x));

        od;

    od;

    return total;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSumsOfDifferenceSetPair( <G>, <N1>, <N2>, <difsums>, <k1>, <k2>, <e>, <lambda> )
##
# FEATURE SomeRefinedDifferenceSumsOfDifferenceSetPair(G, N1, N2, difsums, k1, k2, e, lambda)
InstallGlobalFunction( SomeRefinedDifferenceSumsOfDifferenceSetPair, function (G, N1, N2, difsums, k1, k2, e, lambda)
    local v, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts1, opts2, opt1, opt2, perms1, perms2, perm1, perm2, s1, s2, min, mins, idx;
    if Length(difsums) = 0 then return []; fi;
    v := Size(G);
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);

    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    id := lambda*Size(N2) + e - lambda;
    nonid := lambda*Size(N2);

    newdifsums := []; # will store found difsums

    # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>|G|=",Size(G), ",|N1|=", Size(N1), ",|N2|=", Size(N2), "\n");
    # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>CosetsToCosetsTable(G, N1, N2)=", cosets,"\n");
    for S in difsums do
        # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>S=", S, "\n");
        opts1 := DifferenceSumPreImagesOptions(u3, w2, S[1]);
        opts2 := DifferenceSumPreImagesOptions(u3, w2, S[2]);
        # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>opts1=",opts1,",opts2=",opts2, "\n");
        for opt1 in IteratorOfCartesianProduct(opts1) do
            for opt2 in IteratorOfCartesianProduct(opts2) do
                # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>opt1=",opt1,",opt2=",opt2, "\n");
                # testing identity coeff can be done before permutations
                if not Flat(opt1) * Flat(opt2) = id then
                    continue;
                fi;

                # iterate through all permutations of this partition
                perms1 := DifferenceSumPreImagesPermutationsForced(opt1);
                perms2 := DifferenceSumPreImagesPermutationsForced(opt2);
                # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>perms1=",perms1,",perms2=",perms2, "\n");

                for perm1 in IteratorOfCartesianProduct(perms1) do
                    s1 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm1);
                    for perm2 in IteratorOfCartesianProduct(perms2) do
                        s2 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm2);

                        # continue if s1[idx] <= s2[idx] when k1 > k2 = e
                        if k1 > k2 and k2 = e and s1 <= s2 then continue; fi;

                        # continue when k1 = k2 = e and s1 <> s2
                        if k1 = k2 and k2 = e and s1<>s2 then continue; fi;

                        # forcing choose from the first coset
                        if e > 0 and (s1[1] = 0 or s2[1] = 0) then continue; fi;
                        if e > 1 and (s1[1] < 2 or s2[1] < 2) then continue; fi;

                        # additional check when e > 0
                        if e > 0 then
                            mins := 0;
                            for idx in [1..u2] do
                                min := s1[idx];
                                if min > s2[idx] then
                                    min := s2[idx];
                                fi;
                                mins := mins + min;
                            od;

                            if mins < e then
                                continue;
                            fi;
                        fi;

                        # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>perm1=",perm1,",perm2=",perm2, "\n");
                        # Print("SomeRefinedDifferenceSumsOfDifferenceSetPair>>>s1=",s1,",s2=",s2, "\n");
                        if IsDifferenceSumByTableOfDifferenceSetPair(dif, s1, s2, u2, nonid) then
                            Add(newdifsums, [s1,s2]);
                            # TODO break?
                            break;
                        fi;
                    od;

                od;
            od;
        od;
    od;
    return newdifsums;
end );


#############################################################################
##
#F  SomeRefinedDifferenceSumsOfDifferenceSetPairOptimized( <G>, <N1>, <N2>, <difsums>, <k1>, <k2>, <e>, <lambda> )
##
# FEATURE SomeRefinedDifferenceSumsOfDifferenceSetPairOptimized(G, N1, N2, difsums, k1, k2, e, lambda)
InstallGlobalFunction( SomeRefinedDifferenceSumsOfDifferenceSetPairOptimized, function (G, N1, N2, difsums, k1, k2, e, lambda)
    local v, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts1, opts2, opt1, opt2, perms1, perms2, perm1, perm2, s1, s2, min, mins, idx, multiplier, tracksOfGroup, track, flag, parts1, parts2;
    if Length(difsums) = 0 then return []; fi;
    v := Size(G);
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);

    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    id := lambda*Size(N2) + e - lambda;
    nonid := lambda*Size(N2);

    newdifsums := []; # will store found difsums

    # optimization：在同一轨道的陪集，其差集和元素相同
    multiplier := Multipliers(G, e, lambda)[1];
    tracksOfGroup := TracksOfGroup(G/N2, multiplier);

    for S in difsums do
        opts1 := DifferenceSumPreImagesOptions(u3, w2, S[1]);
        opts2 := DifferenceSumPreImagesOptions(u3, w2, S[2]);
        for opt1 in IteratorOfCartesianProduct(opts1) do
            for opt2 in IteratorOfCartesianProduct(opts2) do
                if not Flat(opt1) * Flat(opt2) = id then
                    continue;
                fi;

                perms1 := DifferenceSumPreImagesPermutationsForced(opt1);
                perms2 := DifferenceSumPreImagesPermutationsForced(opt2);

                for perm1 in IteratorOfCartesianProduct(perms1) do
                    s1 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm1);
                    for perm2 in IteratorOfCartesianProduct(perms2) do
                        s2 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm2);

                        # continue if s1[idx] <= s2[idx] when k1 > k2 = e
                        if k1 > k2 and k2 = e and s1 <= s2 then continue; fi;

                        # continue when k1 = k2 = e and s1 <> s2
                        if k1 = k2 and k2 = e and s1<>s2 then continue; fi;

                        # forcing choose from the first coset
                        if e > 0 and (s1[1] = 0 or s2[1] = 0) then continue; fi;
                        if e > 1 and (s1[1] < 2 or s2[1] < 2) then continue; fi;

                        # additional check when e > 0
                        if e > 0 then
                            mins := 0;
                            for idx in [1..u2] do
                                min := s1[idx];
                                if min > s2[idx] then
                                    min := s2[idx];
                                fi;
                                mins := mins + min;
                            od;

                            if mins < e then
                                continue;
                            fi;
                        fi;

                        flag := true;
                        for track in tracksOfGroup do
                            if Length(track)=1 then
                                continue;
                            else
                                parts1 := s1{track};  # 该轨道下对应的差集和元素，理论上应该相等
                                parts2 := s2{track};
                                if(Length(Set(parts1))>1 or Length(Set(parts2))>1) then
                                    flag := false;
                                    break;
                                fi;
                            fi;
                        od;
                        if flag = true and IsDifferenceSumByTableOfDifferenceSetPair(dif, s1, s2, u2, nonid) then
                            Add(newdifsums, [s1,s2]);
                            break;
                        fi;
                    od;

                od;
            od;
        od;
    od;
    return newdifsums;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSumsOfDivisibleDifferenceSet( <G>, <N1>, <N2>, <difsums>, <w>, <v>, <l>, <lambda1>, <lambda2> )
##
# FEATURE SomeRefinedDifferenceSumsOfDivisibleDifferenceSet(G, N1, N2, difsums, w, v, l, lambda1, lambda2)
InstallGlobalFunction( SomeRefinedDifferenceSumsOfDivisibleDifferenceSet, function (G, N1, N2, difsums, w, v, l, lambda1, lambda2)
    local n, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts, opt, perms, perm, s;
    if Length(difsums) = 0 then return []; fi;
    n := w * v;
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);
    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    # id := l + lambda1*(v-1) + lambda2*(n-v);
    # nonid := lambda1*v + lambda2*n;  # n = l - lambda1 - v * lambda2
    id := l - lambda1 + lambda1 * (v / Size(N2)) + lambda2 * Size(N2) - lambda2 * (v / Size(N2));
    nonid := lambda1 *(v / Size(N2)) + lambda2 * Size(N2);

    newdifsums := []; # will store found difsums

    for S in difsums do
        opts := DifferenceSumPreImagesOptions(u3, w2,  S);
        for opt in IteratorOfCartesianProduct(opts) do
            if not Sum(Flat(opt), x->x^2) = id then
                continue;
            fi;

            # iterate through all permutations of this partition
            perms := DifferenceSumPreImagesPermutationsForced(opt);
            for perm in IteratorOfCartesianProduct(perms) do
                s := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm);
                if IsDifferenceSumByTableOfDivisibleDifferenceSet(dif, s, u2, nonid) then
                    Add(newdifsums, s);
                fi;

            od;

        od;

    od;

    return newdifsums;
end );

#############################################################################
##
#F  SomeRefinedDifferenceSumsOfDivisibleDifferenceSetPair( <G>, <N1>, <N2>, <difsums>, <w>, <v>, <l1>, <l2>, <le>, <lambda1>, <lambda2> )
##
# FEATURE SomeRefinedDifferenceSumsOfDivisibleDifferenceSetPair(G, N1, N2, difsums, w, v, l1, l2, le, lambda1, lambda2)
InstallGlobalFunction( SomeRefinedDifferenceSumsOfDivisibleDifferenceSetPair, function (G, N1, N2, difsums, w, v, l1, l2, le, lambda1, lambda2)
    local n, dif, cosets, w2, u1, u2, u3, id, nonid, newdifsums, S, opts1, opts2, opt1, opt2, perms1, perms2, perm1, perm2, s1, s2;
    if Length(difsums) = 0 then return []; fi;
    n := w*v;
    dif := DifferenceTable(G/N2);
    cosets := CosetsToCosetsTable(G, N1, N2);
    w2 := Size(N2);
    u1 := Size(G/N1);
    u2 := Size(G/N2);
    u3 := Size(N1/N2);
    # id := le + lambda1*(v-1) + lambda2*(n-v);
    # nonid := lambda1*v + lambda2*n;  # n = le - lambda1 - v * lambda2
    id := le - lambda1 + lambda1 * (v / Size(N2)) + lambda2 * Size(N2) - lambda2 * (v / Size(N2));
    nonid := lambda1 *(v / Size(N2)) + lambda2 * Size(N2);

    newdifsums := []; # will store found difsums

    for S in difsums do
        opts1 := DifferenceSumPreImagesOptions(u3, w2,  S[1]);
        opts2 := DifferenceSumPreImagesOptions(u3, w2,  S[2]);
        for opt1 in IteratorOfCartesianProduct(opts1) do
            for opt2 in IteratorOfCartesianProduct(opts2) do
                if not Flat(opt1) * Flat(opt2) = id then
                    continue;
                fi;
                perms1 := DifferenceSumPreImagesPermutationsForced(opt1);
                perms2 := DifferenceSumPreImagesPermutationsForced(opt2);
                for perm1 in IteratorOfCartesianProduct(perms1) do
                    s1 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm1);
                    for perm2 in IteratorOfCartesianProduct(perms2) do
                        s2 := DifferenceSumPreImagesTranslate(cosets, u1, u2, u3, perm2);
                        if IsDifferenceSumByTableOfDivisibleDifferenceSetPair(dif, s1, s2, u2, nonid) then
                            Add(newdifsums, [s1,s2]);
                        fi;
                    od;

                od;
            od;
        od;
    od;

    return newdifsums;
end );

#############################################################################
##
#F  DifferenceSumPreImagesOptions( <u>, <w>, <S> )
##
InstallGlobalFunction( DifferenceSumPreImagesOptions, function (u, w, S)
    local opts, i, parts;

    opts := [1..Length(S)];
    for i in [1..Length(S)] do
        # partition the large coset coeff into the new small coset coeffs
        # adjust call to RestrictedPartitions since it won't use 0 as a summand
        parts := RestrictedPartitions(S[i]+u, [1..w+1], u);
        parts := List(parts, x->List(x, y->y-1));
        opts[i] := parts;
    od;

    return opts;
end );

#############################################################################
##
#F  DifferenceSumPreImagesPermutations( <opt> )
##
InstallGlobalFunction( DifferenceSumPreImagesPermutations, function (opt)
    local perms;

    perms := List(opt, x->PermutationsList(x));

    return perms;
end );

#############################################################################
##
#F  DifferenceSumPreImagesPermutationsForced( <opt> )
##
InstallGlobalFunction( DifferenceSumPreImagesPermutationsForced, function (opt)
    local perms;

    perms := List(opt, x->PermutationsList(x));

    # if the trivial coset contains anything we can force pick 1
    if Sum(perms[1][1]) > 0 then
        perms[1] := Filtered(perms[1], x->(x[1] > 0));
    fi;

    return perms;
end );

#############################################################################
##
#F  DifferenceSumPreImagesTranslate( <cosets>, <u1>, <u2>, <u3>, <perm> )
##
InstallGlobalFunction( DifferenceSumPreImagesTranslate, function (cosets, u1, u2, u3, perm)
    local S, i, j;

    # map the small coset values from the permutation to the appropriate
    # places in the list representing the new potential difference sum
    S := List([1..u2]);
    for i in [1..u1] do
    for j in [1..u3] do
        S[cosets[i][j]] := perm[i][j];
    od;
    od;

    return S;
end );

#############################################################################
##
#F  IsDifferenceSumByTable( <table>, <S>, <v>, <k>, <lambda>, <w> )
##
InstallGlobalFunction( IsDifferenceSumByTable, function (table, S, u, nonid)
    # method is identical to IsDifferenceSum but uses precomputed values
    local i, j, index, count;

    count := List([1..u], x->0);

    for i in [1..u] do
    for j in [1..u] do

        index := table[i][j];
        count[index] := count[index] + S[i]*S[j];

        if not index = 1 and count[index] > nonid then
            return false;
        fi;

    od;
    od;

    return true;
end );

#############################################################################
##
#F  IsDifferenceSumByTableOfDifferenceSetPair( <table>, <S1>, <S2>, <u>, <nonid>)
##
# FEATURE IsDifferenceSumByTableOfDifferenceSetPair(table, S1, S2, u, nonid)
InstallGlobalFunction( IsDifferenceSumByTableOfDifferenceSetPair, function (table, S1, S2, u, nonid)
    local i, j, index, count;

    count := List([1..u], x->0);

    for i in [1..u] do
    for j in [1..u] do

        index := table[i][j];
        count[index] := count[index] + S1[i]*S2[j];

        if not index = 1 and count[index] > nonid then
            return false;
        fi;

    od;
    od;

    return true;
end );

#############################################################################
##
#F  IsDifferenceSumByTableOfDivisibleDifferenceSet( <table>, <S>, <u>, <nonid> )
##
# FEATURE IsDifferenceSumByTableOfDivisibleDifferenceSet(table, S, u, nonid)
InstallGlobalFunction( IsDifferenceSumByTableOfDivisibleDifferenceSet, function (table, S, u, nonid)
    local i, j, index, count;

    count := List([1..u], x->0);

    for i in [1..u] do
    for j in [1..u] do

        index := table[i][j];
        count[index] := count[index] + S[i]*S[j];

        if not index = 1 and count[index] > nonid then
            return false;
        fi;

    od;
    od;

    return true;
end );

#############################################################################
##
#F  IsDifferenceSumByTableOfDivisibleDifferenceSetPair( <table>, <S1>, <S2>, <u>, <nonid> )
##
# FEATURE IsDifferenceSumByTableOfDivisibleDifferenceSetPair(table, S1, S2, u, nonid)
InstallGlobalFunction( IsDifferenceSumByTableOfDivisibleDifferenceSetPair, function (table, S1, S2, u, nonid)
    local i, j, index, count;

    count := List([1..u], x->0);

    for i in [1..u] do
    for j in [1..u] do

        index := table[i][j];
        count[index] := count[index] + S1[i]*S2[j];

        if not index = 1 and count[index] > nonid then
            return false;
        fi;

    od;
    od;

    return true;
end );

#############################################################################
##
#E

